arXiv: 1507.03448v6 [math.NA] 3 Aug 2016 


Augmenting Numerical Stability of the Galerkin Finite Element Formulation 
for Electromagnetic Flowmeter Analysis 

Sethupathy S. 1 *, Udaya Kumar 2 

1,2 Department of Electrical Engineering, Indian Institute of Science, Bangalore 560012, India 
*sethupathys87@ gmail.com 


This paper is a preprint of a paper accepted by 1ET Science, Measurement & Technology and is subject 
to Institution of Engineering and Technology Copyright. When the final version is published, the copy of 
record will be available at IET Digital Library. 



Abstract: Due to the complexities in handling liquid metals, theoretical evaluation of the sensitivity 
of magnetic flow meters forms an attractive and preferred choice. The classical Galerkin finite 
element formulation is generally opted for the required evaluation. However, it is known to lead to 
numerical oscillations at higher flow rates. To overcome this, modified methods like upwind/Petrov- 
Galerkin schemes are generally suggested in allied areas like fluid dynamics. However, it requires 
the evaluation of stabilization parameter and this parameter is not readily available for elements of 
order beyond quadratic. After a careful analysis of the numerical instability through a reduced one 
dimensional problem, an elegant and stable approach is devised. In this scheme, the input magnetic 
field is restated in terms of the associated vector potential and the classical Galerkin Finite Element 
Method (GFEM) is employed without any modification. The analytical solution of the associated 
difference equation is employed to show: (i) the stability of the proposed approach at higher flow 
rates and (ii) quantification of the small oscillations what it introduces at intermediate flow rates. 
It is then applied to the original flowmeter problem and the stability of the numerical solution is 
clearly demonstrated. 


1. Introduction 

Electromagnetic flowmeter is a non-invasive instrument which is widely used in fast-breeder reactors for 
the measurement of flow rate of liquid metals. As an accurate measurement of flow rate is essential for the 
safe operation and control of the reactor, the performance of flowmeter needs to be reliably ascertained. Due 
to the practical difficulties in handling liquid metals, experimental determination of flowmeter sensitivity is 
an involved job. Hence, accurate theoretical or numerical evaluation of the sensitivity formed an attractive 
alternative. 

Fig-D shows the schematic of electromagnetic flowmeter. The measurement probes (Vj & Vo) are 
placed perpendicular to both magnetic field and flow direction. Circulating currents (Ji, J 2 , J3) are due 
to the spatial variation in the induced electric field. The reaction magnetic field (b rc ) produced by these 
currents cancels the applied magnetic field (B ap ) in the upstream region and aids it in the downstream 
side. This cross-magnetizing effect apparently shifts the effective magnetic field along the flow direction. 
Generally in the analysis, the applied or the ambient and the reaction magnetic fields are separated. 

The governing equations in terms of the magnetic vector potential A of the reaction field and the 
electric scalar potential 0 arising out of current flow are given by [1], [j2]]: 

V ■ (<rV0) — V-(auxVxA) = V-(crux B ap ) (1) 

aV0 — —V 2 A -uuxVxA = (tux B ap (2) 

where, /1 is the magnetic permeability, a is the electrical conductivity and u is the velocity function of the 
fluid flow. The relative strength of the reaction magnetic field (b rc ) to the applied magnetic field (B ap ) is 
indicated by the magnetic Reynolds number R m = fiau z Dh where, Dh is the hydraulic diameter of the 
pipe 0, @. 

For problems with R rn < 1, the reaction magnetic field (b rc ) and hence the cross magnetizing effect 
becomes negligible. In such a situation, the last term on the Left Hand Side (LHS) of (Tj) becomes 


negligible and the resulting equation can be solved independently. It has been reported in [1] that whenever 
the length of the magnet is > 1.5 times the pipe diameter, a two dimensional (2D) analysis across the 
cross section perpendicular to the flow is permissible. For such cases the analytical solution can be found 
[3], [4[. In liquid metals, however, the conductivity is very high and hence the induced currents are large, 
which leads to strong cross-magnetizing effects. As a result, a full three dimensional (3D) analysis will be 
required. Due to the complexity in handling reaction field and the flow-geometry, numerical techniques 
like Galerkin Finite Element Method (GFEM) is generally employed. 



Fig. 1. Schematic of Electromagnetic flowmeter. j3J 


It is well known that GFEM, similar to the central-difference scheme, is averaging in nature and hence it 
is diffusive (6].It is shown to become numerically unstable whenever the convection starts dominating over 
the diffusion. This numerical instability is widely addressed in the fluid dynamics literature dealing with 
transport equation [6j. Streamline Upwind Petrov Galerkin (SUPG) scheme |7], Galerkin Least Squares 
(GLS) |8J, Finite Increment Calculus (FIC) |5] and Multiscale scheme [ lOi] are suggested for stabilizing 
the solution. The same upwinding schemes have also been adopted for electromagnetic problems. For 


example the moving conductor problem has been analyzed in |1T)-[151 using the upwinding Petrov- 
Galerkin scheme. 

A similar numerical instability is also encountered with the GFEM simulation of the magnetic flowmeter. 
The SUPG scheme has been successfully employed in situations involving magnetic Reynolds number 
in the range 4 to 28 0. However, its performance for higher magnetic Reynolds number is yet to 
be quantified. Apart from this, SUPG involves more computation for higher order elements and also 


it is difficult to find the stabilization parameters for elements with order higher than quadratic [16), 
ED- The present work basically aims to overcome these difficulties in the finite element simulation of 
electromagnetic flowmeters. 

The investigations on the numerical instability, as well as, that for the remedial measures, were always 
carried out with the one dimensional version of the original problem j6], [ 161, [18]. This is because, the 
required analysis in 2D or 3D is almost impractical and further, the one dimensional (ID) version of the 
fluid dynamics/thermal problems contained all the required features of the original problem. 

Following the same philosophy, the required analysis for the present work will be carried out using ID 













version of the problem. Both finite difference and Z-transform approaches will be employed for the analysis 
of numerical instability. From an insight thus obtained, a novel and stable scheme will be proposed for 
the classical GFEM. Subsequently the proposed method will be applied to the original flowmeter problem 
and the stability of the scheme will be demonstrated. 


2. Present work 


For the theoretical investigation on the source of numerical instability, the analytical solution of the 
global set of equations given by GFEM is required. However, such an exercise is nearly impractical to 
be carried out in 2D or 3D and hence it is customary to resort to a ID version of the problem 0, [ 161, 

0 . 


Following the same, a reduced one dimensional problem is considered. For this, conducting fluid is 
assumed to occupy the whole space with a spatially uniform velocity u~ in the ^-direction. The magnetic 
field is applied only in x-direction and is defined by, 


{ 0 0 < z < a 

B a < z <b 

0 b < z < L 


(3) 


where, L is the length of the analysis domain and B x exists between a and b. 

With these imposed conditions, the field variables cease to have any variation along x and y directions. 
The governing equations (TJ) & ([2]) therefore reduces to, a single ordinary differential equation (ODE) in 
terms of A y , the only non-vanishing component of A: 


d 2 A 


T y>(ju z 


dA.„ 


f-U7Xl z B x 


(4) 


dz 2 dz 

Left hand side of the equation 0 has the same structure as the one used in fluid dynamics literature 
for investigating the numerical instability [6j, |7|, |9j. 


A. Analysis on instability 

Application of the Galerkin finite element method (GFEM/FEM) and the finite difference method (FDM) 
to 0, both results in same set of difference equation 0, [jl9:|. 

The difference equation for n th node, 


( — 1 — Pe)A y (n-\) + 2 Ay^ + (—1 + Pe)A y (n + i) — 2PeAz B x ^ n ) 


(5) 


where the Peclet number, Pe = (/j,ou z Az) /2 and Az represents the element length. It may be worth 
noting here that, the oscillatory property of the numerical solution is basically a function of Pe rather than 
its constituent parameters taken in isolation. Hence the allied literature considers Pe as the independent 
variable in the required study 0, [H5J, [17]. When Pe > 1, a root of the above difference equation 
becomes negative and it has been identified to be the source of numerical oscillation [20], ED- 

The instability problem can also be analyzed by bringing the tools from the control system theory. The 
difference equation is transformed to the frequency domain for an easier analysis. 







The z-transform of ([5]), 


— 1 — Pe)Z 1 + 2 + (—1 + Pe)Z j A y = 2PeAzB x 


can be written in transfer function form, 
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( 6 ) 


(7) 


B x —1 + Pe (Z — 1) [Z + 1) 

The above transfer function has poles at -1 & +1. The pole located at ’-1’ is responsible for numerical 


oscillations [211. This observation is not specific to any particular excitation. 

In control systems, the controller design is always coupled with the pole-zero cancellation. Zeros of the 
transfer function ([7]) can be modified by changing the input to a suitable form. A novel way to bring the 
desired zeros of the transfer function is to express the input magnetic field B x in terms of the associated 
magnetic vector potential A sy . 

Substituting B x = —dA sy /dz in the equation (|4]), in the difference form: 


( 1 + 2 Ay^ n ) T ( 1 T Pa) A y ^ n _Pe(^A S y^ n _i'j A sy ^ n _ |_i)) 


( 8 ) 


The corresponding transfer function is: 


when Pe » 1 
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( 10 ) 


( 11 ) 


A perfect pole-zero cancellation occurs for Pe » 1, which ensures absolute stability. For very low 
Peclet numbers (< 1), the intrinsic accuracy of GFEM is left unaltered. However, for Peclet numbers in 
the range 1 — 30, the pole-zero cancellation is not perfect and hence some oscillation can arise. In order 
to find the maximum amplitude of this oscillation and the corresponding value of Pe, further work is 
carried out. 


B. Boundary conditions for the one dimensional version of the problem 

The analytical and numerical solutions of difference equations ([5]) & ([8]) will be considered. The input 
parameters employed are p = 47t x 10~ 7 Hm~ l , a = 7.21 x 10 6 Srrr 1 , u z and discretisation length (Az) 
















are kept as variables to make a study on the oscillations. The applied magnetic field B x is presented in 
fig. [2a] along with the corresponding vector potential A sy . 

Because of the infinite spatial extension of the current in the one dimensional version of the problem, 
it cannot be directly mapped on to the original problem. Therefore, it becomes necessary to arrive at the 
appropriate boundary conditions for the same. 

For this purpose, a two dimensional version of the problem is considered, which deemed to have 
adequate representation of the original problem. As shown in the fig. [2bJ except for the fact that the fluid 
is now confined to the region between two finitely conducting plates kept parallel to the xz plane, all 
other aspects are similar to ID problem. The velocity of the fluid is assumed to be uniform in the gap. 
In order to apply the correct boundary conditions for the reaction magnetic field, the air region outside 
the plates up to a distance of 5 d in y direction (where, d is the separation distance between the plates) 
is included in the analysis. The axial length of the computational domain is set to lOd For the FEM 
discretisation linear quadrilateral elements are employed. The element size along the axial direction is 
kept constant while that along the y direction is varied to accurately capture the current flow near the 
pipe wall. The element size selected always satisfied Pe < 1. 

Even though the spatial extension of the current in the 2D problem is also infinite in the x direction, 
the induced currents form closed loops in the yz plane. Hence, the boundary condition (A = 0) which is 
relevant to the original problem is enforced at far distances of the upstream and downstream sides (fig. 


2b). For a quantitative assessment of the reaction magnetic field and induced currents, numerical solution 


of the 2D problem is sought. 

Selected results from the simulation are presented in fig. [2c] & 2d for different gaps (d) between the 
plates. In line with the expectation, it was found that for a wide range of velocities and the separation 
gap between the plates, the total current crossing the bisecting plane (shown in fig. ( [2b] )) is zero. This 
ensured that the reaction magnetic field b x vanishes at far distances along the downstream side. Further 
the increase in the gap ( d ) between the plates, b x is almost confined close to the region where the applied 
magnetic field exists. 

The solution for the 2D version of the problem has clearly shown that magnetic field is dragged only 
along the downstream side, while it is always confined at the upstream end (for the range of velocities of 
interest, which correspond to yau z > 10). Based on this for the ID problem, A y is specified to be zero 
at the upstream boundary and the vanishing b x leads to dA y /dz = 0 at the downstream boundary. 

Along with these boundary conditions and the input magnetic field defined in ([3]), the analytical solution 
of the ODE ([4]) is: 
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( 12 ) 


where, k = yau z \ B = value of B x in a < z < b as mentioned in ([3]). This can serve as the reference 








(C) (d) 

Fig. 2. (a) Input quantity and its derivative form, (b) Schematic of 2D problem, (c) 2D: Computed reaction magnetic field for different plate 
separation distances ( d ) with u z = 10 ms _1 . (d) 2D: Computed magnetic vector potential for the same. 


for the evaluation of the error in the numerical solution of the governing equation. Sample results obtained 
from the FEM are presented in fig. [3] along with the analytical solution of the governing equation, where 
the reaction magnetic field b x is calculated using the forward difference scheme: 

Ay(n + 1) - A y (n) 
b x (n) = -—- i - -^- L 

The solution is stable for Pe » 1 (practically verified for Pe ranging from 30 to 30000) with A sy as 
input. However, as mentioned earlier, in the mid-range 1 < Pe < 30 oscillation exists. Sample results are 
presented in fig. [3] 

In order to identify the location (the value of Pe) and peak amplitude of the oscillation, analytical 
solution of the FEM equation, which appears in the difference equation form, is performed in the next 
section. 

C. Location and value of the maximum error 

The governing difference equations can be rewritten as: 

For magnetic field B x as the input: 


rA y (n — 1 ) + (—1 — r)A y (n) + A y (n + 1 ) = 2PeAz B x (n ) 


( 13 ) 
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Fig. 3. FDM/FEM solution, (a) Pe = 3000, A« = 0.25. (b) Pe = 3, A* = 0.25. 


Similarly for vector potential A sy as the input 


rAy(n - 1) + (-1 - r)A y (n) + Mn + 0 = - 1) - A„(n + 1)) (14) 

where, r = (—1 — Pe )/(—1 + Pe) ; 

Using the method described in [20], the general solution of the equations (fl~3|) and ([14]) can be found 


as, 


A y (n) = + k 2 r n + y p (n ) 


(15) 


where, y p (n) is the particular solution and k \, k 2 are parameters of the complementary solution. 

The direct solution of the difference equation with piecewise input is difficult and hence the problem 


domain is divided into five sub-domains ( B , C, D, F & G) as shown in figs. 7a & 8a of Appendix. The 
solution steps are also detailed in the Appendix along with its validation. 

The oscillation in the numerical result can be found both analytically and numerically to occur at the 
end of domain C. Based on this, for finding the peak amplitude of the oscillation and corresponding Pe, 


solution for domain C i.e., equation (311 is considered. The reaction magnetic field at the end of domain 

A y (m c — 1) — A y (m c ) c 2 (r mc 1 — r mc ) A 

Az 


b n = 


Az 


Az 


(16) 


Equation (fl6|), has a constant part given by A Az and an oscillatory part given by (c 2 (r mc 1 — r mc ))/Az. 








































Comparing with the analytical solution given in equation ( fl2| ), oscillatory part can be identified to be 
the major part of the error in the numerical solution. For the range of velocities of interest (i.e. jiau z > 10) 
this part amounts to almost the total error. 


Let, 



^ C2 (j-m c — 1 _ r m c \ 

ba ~ Az 

(17) 

substituting (50) in (|17|), 

r 1 + r 5(1 - Pe) 

ba 2r 2 ~ (1 + Pe) 2 

(18) 

Similarly, for B x input, 

r B 5(1 - Pe) 

r 1 + Pe 

(19) 



Pe 


Fig. 4. % peak error in the numerical solution for a wide range of Peclet numbers. 


The magnitude of the oscillation present in the analytical solution of the difference equation is given 
by (18) & ( p~9] ) for the respective inputs (i.e., A sy & B x respectively). The error obtained from the above 
is plotted in fig. |4j The peak error, with input field specified in terms of the vector potential occurs at 
Pe = 3 and its magnitude in % is 1/8. 

The above numerical exercise has once again confirmed that the proposed scheme is very stable for 
large flow rates. Also, even in the midrange of flowrates (and hence Pe), the error in the numerical results 
is lower than that with GFEM and further it can be controlled by opting for different discretisation. 


D. Performance with quadratic elements 

Unlike that with the first order elements, SUPG scheme for higher order elements requires more 
computation. On the other hand, the proposed scheme is free of such issues and when implemented 
for second order elements, its performance is found to be equally good. Sample numerical results are 
presented in fig. [5] as an evidence for the same. 
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Fig. 5. Verification for quadratic elements, (a) Pe = 3000. (b) Pe = 3. 
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E. Analysis for the flowmeter 

Up till now all the analysis was limited to one dimensional version of the problem and therefore it 
was deemed necessary to scrutinize the proposed scheme with the original problem. For this, governing 
equations 0 & (|2j) are solved using GFEM [|2j. The input magnetic field B ap is set to have only 
the ^-component B x , whose magnitude varies only along the flow direction (as shown in fig [6a]). The 
corresponding vector potential A s will have only y component. It is readily available if one employs 
numerical simulation for the input magnetic field or it can be obtained from —dA sy /dz = B x when the 
measured flux density is employed. 

For the numerical experiment, a non-magnetic stainless steel pipe with inner and outer diameters 
517 mm & 557 mm respectively is considered. The conductivity of the pipe material is 1.16 x 10 6 Sm~ l 
and that of the liquid sodium inside is 7.21 x 10 6 Sm~ l . It is true that due to cavitation and other associated 
problems, velocities beyond few to few tens of meters are impractical with liquid metals. Nevertheless, 
simulations are carried out for velocities ranging up to 3000 ms~ l (which corresponds to R m = 14851), 
solely to demonstrate the robustness of the proposed scheme. The intention was to consider the possible 
application of the proposed scheme for general moving conductor problems. Incidentally, R m is related 
to the Peclet number through Pe = R m Az/2D h . 

Sample simulation results are compared with that obtained from the original formulation in fig[6] The 
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as inputs, (c) x-component of reaction magnetic field b r 
b rc in x = 0 plane with V x A s as input. 


x = 0 plane with B ap as input. 


& V x A s 

(d) ^-component of reaction magnetic field 


reaction magnetic field for Pe = 40 as given by two approaches is presented in fig. 6c & 6d while 


fig. 6b presents reaction magnetic field along the pipe axis for different values of Pe. From the results 
obtained from extensive simulations, it is observed that for the midrange of Pe there exists a small amount 
of oscillation, which asymptotically vanishes with increase in Pe. These observations are in line with the 
inference drawn earlier from the ID analysis. It must be noted here that all the analysis reported here 
considers only an axially varying applied magnetic field, which generally holds true for flowmeters. 


3. Summary and Conclusion 

Theoretical evaluation of the sensitivity of electromagnetic flowmeter seems to be the best possible 
choice especially when it is to be used for the measurement of liquid metal flows. The commonly employed 
Galerkin finite element formulation is known to become unstable for larger flow rates. SUPG scheme 
is generally suggested in the pertinent literature for overcoming this problem. However SUPG scheme 
requires computation of the stabilization parameter, which involves a lot more calculations with higher 
order elements. In addition, it is difficult to arrive at the stabilization parameters for elements with order 
beyond quadratic [ [T6| , |T7J. 




















































By analyzing the one dimensional version of the problem a novel scheme has been devised which is 
free of above mentioned difficulties. In this scheme, classical GFEM is retained intact and only the input 
magnetic field is restated in terms of the associated vector potential. The analytical solution of the FEM, 
for the proposed scheme indicate that solution exhibits a small amount of oscillation in the mid-range 
of flow rates (1 < Pe < 30) and this oscillation asymptotically vanishes with the increase in Pe (for 
Pe > 3). The maximum error due to this numerical oscillation has been quantified analytically. Even the 
maximum error, which is shown to occur at Pe = 3, is lower than what is found with flux density as 
the input. Finally, the proposed scheme is applied to the original flowmeter problem and the simulation 
results indicate that inferences drawn on the ID version of the problem remain valid even for the 3D 
case. In summary, a simple and robust scheme has been proposed for the FEM solution of the flowmeter 
problem involving only an axially varying applied magnetic field. 
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Fig. 7. (a) Five sub-domains and their ranges — B x as Input, (b) Validation of the analytical solution with B x as input Pe = 200, Az = 
0.25, m c = 6, = 31, = 31. (c) Validation of the analytical solution with B x as input Pe = 3, Az = 0.33, m c = 4, trtt = 

23, m d = 23. 
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Appendix 

A. Analytical solution - B x as Input 

The direct solution of a difference equation with piecewise-defined input is difficult. Hence the domain 
is divided into five sub-domains as shown in fig. [7a] 

Following notations A = B A z, r = (—1 — Pe)/(— 1 + Pe), and y = A y are employed in the subsequent 
steps. In order to distinguish the solutions of each domain, the parameters of complementary solution are 
designated with their domain names (fq, b 2 , Cy, c 2 etc.,). And the other variables such as y, particular 
solution y p and node number n are suffixed with their corresponding domain names (y b , rib. He, n c , y pc 
etc., ). The governing difference equations and their general solutions for different domains are listed 
below. 

For domain B (0 < iib < m b ): 


ry b (n b - 1) + (-1 - r)y b [n b ) + y b (n b + 1) = 0 


(20) 




























=>■ Vb{n b ) = 6 i + 6 2 r nb 


( 21 ) 


For domain F (0 < rij < 1): 

ry f (n f - 1) + (-1 - r)y f (n f ) + y f {n f + 1) = A(1 - r)n f (22) 

=> VfM = fi + / 2 r n f + y vf (yif) (23) 

For domain C (0 < n c < m c ): 

ry c {n c - 1) + (-1 - r)y c (n c ) + y c {n + 1) = A(1 - r) (24) 

=>■ y c (n c ) = Cl + c 2 r nc + y pc {n c ) (25) 

For domain G (0 < n g < 1): 

ry g (n g - 1) + (-1 - r)y g (n g ) + y g (n g + 1) = \{l-r)(l-n g ) (26) 

=> y g {n g ) = gi + g 2 r ng + y pg (n g ) (27) 

For domain D (0 < < m^): 

ryd(nd - 1) + (-1 - r)y d (n d ) + y d (n d + 1) = 0 (28) 

=> Vd(n d ) = rfi + d 2 r nd (29) 


The difference equations in domains F, 67 & G are inhomogeneous and their particular solutions are 


computed [20]. The complete general solutions of domains F, C & G are, 


VfM = fi + hr ns + + \ n ) ( 3 °) 

y c (n c ) = ci + c 2 r nc + A n c (31) 

y 9 ( n g) =91 + 92T na + (1 - 2 (r + -\) ) Xn<} ~ ^ (32) 

The global boundary conditions imposed on the extreme boundaries. It helps in eliminating two 
variables. 

2/ b (0) = 0 =► 61 = - 6 2 (33) 

yd{m d ) = y d {m d + 1) => d 2 = 0 (34) 


The intermittent boundaries should satisfy the equality condition and the governing equation itself at 
the joining nodes. At junction G — D, 

Equality condition: 

A 

= d 1 


1 — r 


2/ ff (l) = 2/d(0) =>• gi + g2 r r + 


(35) 






Satisfying, governing equation 


ry g {0) - (1 + r)y d {0) + y d (l) = 0 gi + g 2 = di (36) 

The above mentioned conditions are employed for other three junctions and six more equations are 


obtained. These six equations are solved along with equations ( [33] ), ( [34] ), ( ]35[ ) & ([36]). The parameters of 
complementary solutions are obtained as: 


9 2 
c 2 
/ 2 

b 2 

h 

Cl 

ffi 

cl 1 


(1 
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r mc (l — r) 

-A- 1 


(1 




A 


+ 


/a 


rpTTlb 1 ^ 

bi + b 2 r m " - f 2 
r 

fi + hr + A-- - c 2 

r — 1 

ci + c 2 r mc + A m c - g 2 
A 

9 1 + W + - - 


(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


The analytical solution of the difference equation is compared with the numerical solution obtained 


from FDM and FEM in fig. 7b & 7c 


B. Analytical solution with derivative of A sy as Input 


The input in terms of A sy is plotted in fig. 8a Similar to B x input, the problem is divided into five 
sub-domains. 

The governing difference equations and their general solutions for domains B , C. & I) are same as 
that of the previous case, with A redefined as A = ( A sy (n c — 1) — A sy (n c + l))/2 and it is equivalent 
to B Az. F & G domain equations are listed as follows, 


For domain F (0 < nj < 2): 


ryf(n f - 1) + (-1 - r)y f [n f ) + y f (n f + 1) = A 


1 — r 


~ n f 


VfM = h + f2T nf + rtf + ^ n) 


For domain G (0 < n q < 2): 


ry g (n g - 1) + (-1 - r)y g (n g ) + y g (n g + 1) = A(2 - n g ) 

^ y g (n g ) =gi + g 2 r n ° + (l - ^ r + _\^ ) Xn a ~ 


1 — r 


(45) 

(46) 

(47) 

(48) 
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Fig. 8. (a) Five sub-domains and their ranges A ay as input, (b) Validation of the analytical solution with A sy as input Pe = 300, Az 
0.33, m c = 4, mb = 22, md = 22. (c) Validation of the analytical solution with A sy as input Pe = 3, Az = 0.25, m c = 6, m;, 
30, md = 30. 


- FDM/FEM simulation 

Analytical solution of the 
difference equation 



The computed values of parameters in the complementary solution are: 
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-A 
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2(r - 1) 
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9 1 + . 92 ?" + 


A(r — 2) 

2(r — 1) 


(49) 

(50) 

(51) 

(52) 

(53) 

(54) 

(55) 

(56) 


The analytical solution of the difference equation is compared with the numerical solution obtained 


from FDM and FEM in fig. 8b & 8c 
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